PD_scRNAseq: Monocle3

Author

Brian M. Schilder, Bioinformatician II

Setup

start_time <- Sys.time() 

#### #### Load Objects & Functions #### #### 
######################################################
# Import functions
root = "./"
source(file.path(root,"MAIN.R"))
import_parameters(params)

## **** __Utilized Cores__ **** = 2$subsetGenes
## [1] "protein_coding"
## 
## $subsetCells
## [1] 500
## 
## $resolution
## [1] 0.6
## 
## $resultsPath
## [1] "./Results"
## 
## $nCores
## [1] 2
## 
## $perplexity
## [1] 30

# load("Results/Current_Pipeline/scRNAseq_results.RData")
load(file.path(resultsPath, "scRNAseq_results.RData")) 
######################################################


#### #### PACKAGES #### #### 
######################################################
print("Written using: Seurat version* 2.3.4 2018-07-17")

## [1] "Written using: Seurat version* 2.3.4 2018-07-17"

# https://satijalab.org/seurat/install.html
library(Seurat) #
paste("Seurat", packageVersion("Seurat")) 

## [1] "Seurat 2.3.4"

# library(monocle) # BiocManager::install("monocle")   
# paste("monocle", packageVersion("monocle")) 
library(monocle3)
paste("monocle3", packageVersion("monocle3")) 

## [1] "monocle3 0.1.2"

library(garnett) # devtools::install_github("cole-trapnell-lab/garnett")
paste("garnett", packageVersion("garnett"))

## [1] "garnett 0.2.3"

library(cowplot) 
library(ggplot2)
library(dplyr) 
library(data.table) 
library(readxl) 
library(reshape2)
library(ggrepel)
library(plotly)
library(GeneOverlap) # BiocManager::install("GeneOverlap") 
######################################################


# Exporting 3D plots 
knitr::knit_hooks$set(webgl = rgl::hook_webgl)
dge_limit <- F # 100
nCores <- 4#parallel::detectCores()
set.seed(2019)

Preprocessing

Seurat/2.3.4 to Monocle3

Also subset to only protein-coding genes.

biotypes <- get_biotypes(DAT)

## Subsetting genes: protein_coding

proteins <- biotypes[biotypes$gene_biotype=="protein_coding",]$hgnc_symbol %>% droplevels()
protDAT <- subset_seurat(DAT, genes.use = proteins)

cds <- seurat_to_monocle(seurat_object = DAT, monocle_version = 3)

## Object representation is consistent with the most current Seurat version.
## [1] "Processing..."
## [1] "+ Expression data"
## [1] "+ Phenotype data"
## [1] "+ Feature Data"
## [1] "Expression Data dims:"
## [1] 24914 19144
## [1] "Metadata dims:"
## [1] 19144    17
## [1] "Feature Data dims:"
## [1] 24914     1
## [1] "+ Converting to monocle (Version 3)"

cds <- cds[levels(proteins),]

Normalization & PCA

cds <- monocle3::preprocess_cds(cds, 
                                num_dim = 20, #100 by default
                                residual_model_formula_str = "~ nUMI + percent.mito")
monocle3::plot_pc_variance_explained(cds) 

Dimensionality Reduction

Using UMAP.

# c("UMAP", "tSNE", "PCA", "LSI")
cds <- monocle3::reduce_dimension(cds, reduction_method = "UMAP", 
                                  max_components = 3
                                  # umap.fast_sgd = T,
                                  # cores=nCores
                                  ) 

Clustering

• Unsupervised clustering of cells is a common step in many single-cell expression workflows. In an experiment containing a mixture of cell types, each cluster might correspond to a different cell type. - This function takes a cell_data_set as input, clusters the cells using Louvain community detection, and returns a cell_data_set with internally stored cluster assignments.

• In addition to clusters this function calculates partitions, which represent superclusters of the Louvain communities that are found using a kNN pruning method. Cluster assignments can be accessed using the clusters function and partition assignments can be accessed using the partitions function.

• Using only the topN variable genes to cluster

Find Variable Genes

# Method 1
# variable.genes <- cds@preprocess_aux@listData$gene_loadings[,1:3] %>%
#   abs() %>% rowSums() %>% base::sort(decreasing = T)
# head(variable.genes) 
varDAT <- Seurat::FindVariableGenes(object = protDAT, 
                             mean.function = ExpMean, 
                             dispersion.function = LogVMR, 
                             selection.method ="dispersion", do.plot = T, 
                             top.genes = 2000)

## Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
## Binning grid too coarse for current (small) bandwidth: consider increasing
## 'gridsize'

var.genes <- varDAT@var.genes
head(var.genes)

## [1] "TMSB4X" "FTL"    "S100A8" "FTH1"   "B2M"    "S100A9"

Louvain Clustering

cds <- monocle3::cluster_cells(cds,
                               resolution = 1e-4,#c(10^seq(-6,-1)), # Do NOT set to high number 
                               reduction_method = "UMAP", 
                               cores = nCores,
                               clustering_genes = var.genes)
# Add cluster info to metadata (for easier DGE analysis)
pData(cds)$Cluster <- monocle3::clusters(cds)
pData(cds)$Partition <- monocle3::partitions(cds)
# 3D plot
p3d <- plot_cells_3d(cds, color_cells_by = "cluster", show_trajectory_graph = F)
p3d

# Save each cluster as a separate OBJ file that can later be 
# imported into a 3D modeling platform (e.g. Blender, Meshlab) for more advanced visualization
## save_3D_clusters(cds)

Pseudotime

# Learn
cds <- monocle3::learn_graph(cds)

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |=================================================================| 100%

## Warning in louvain_clustering(data, pd[row.names(data), ], k = k, weight = weight, : RANN counts the point itself, k must be smaller than
## the total number of points - 1 (all other points) - 1 (itself)!

# Order cells
root_pr_nodes = get_earliest_principal_node(cds, variable = "dx", variable_value = "PD")
cds <- monocle3::order_cells(cds, root_pr_nodes = root_pr_nodes) 

monocle3::plot_cells(cds)

monocle3::plot_genes_in_pseudotime(cds[c("CD14","FCGR3A"), ] )

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

UMAP Plots

Disease status (dx), mutation status (mut), and individual ID (ID) show good mixture across clusters.

for(v in c("dx","mut","ID")){
  cat("\n")
  cat("####",v,"\n")
  p <- monocle3::plot_cells(cds, color_cells_by=v, label_cell_groups=F, show_trajectory_graph = F)
  print(p)
  cat("\n")
} 

dx

mut

ID

Cell-type Identication

Garnett

Use a pre-trained classifier from Pline et al.

library(org.Hs.eg.db)

## Loading required package: AnnotationDbi

## 
## Attaching package: 'AnnotationDbi'

## The following object is masked from 'package:plotly':
## 
##     select

## The following object is masked from 'package:dplyr':
## 
##     select

## 

columns(org.Hs.eg.db)

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [25] "UNIGENE"      "UNIPROT"

# load(url("https://cole-trapnell-lab.github.io/garnett/classifiers/hsPBMC")) 
load("./Data/Garnett/hsPBMC")
cds = garnett::classify_cells(cds, 
                               hsPBMC,
                               db = org.Hs.eg.db,
                               cluster_extend = T,
                               cds_gene_id_type = "SYMBOL")

## Warning in doTryCatch(return(expr), name, parentenv, handler): The
## following genes used in the classifier are not present in the input CDS.
## Interpret with caution. ENSG00000174059The following genes used in the
## classifier are not present in the input CDS. Interpret with caution.
## ENSG00000007062The following genes used in the classifier are not present
## in the input CDS. Interpret with caution. ENSG00000157404The following
## genes used in the classifier are not present in the input CDS. Interpret
## with caution. ENSG00000185291

monocle3::plot_cells(cds,
                     group_cells_by="cluster",
                     color_cells_by="cluster_ext_type",  # cell_type
                     show_trajectory_graph = F )

Save Checkpoint

Save R object and run memory-intensive DGE analyses on computing cluster.

saveRDS(cds, file="./Data/monocle3_CDS.RData")